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Abstract. The Fermi-Pasta-Ulam (FPU) pioneering numerical experiment 
played a major role in the history of computer simulation because it introduced 
this concept for the first time. Moreover, it raised a puzzling question which was 
answered more than 10 years later. After an introduction to this problem, we 
briefly review its history and then suggest some simple numerical experiments, 
with a provided Matlab^ code, to study various aspects of the "FPU" problem. 



1. The physical question 

The "FPU problem", as it is known presently, bears the name of the three scientists || 
who wanted to study the thermalization process of a soHd pp. As revealed by S. 
Ulam later these authors were looking for a theoretical physics problem suitable 
for an investigation with one of the very first computers, the "Maniac". Their work, 
published in 1955 in a classified Los Alamos National Laboratory report, had been 
completed shortly before the death of Enrico Fermi in 1954. Ulam said later that 
Fermi did not suspect the importance of this discovery and considered this work as 
"minor" . He nevertheless intended to talk about it at the conference of the American 
Mathematical Society, where he had been invited just before he became seriously ill. 

Fermi, Pasta and Ulam decided to study how a crystal evolves towards thermal 
equilibrium by simulating a chain of particles of unitary mass, linked by a quadratic 
interaction potential, but also by a weak nonlinear interaction. One of the one- 
dimensional systems they considered is described by the Hamiltonian 

where Ui is the displacement along the chain of atom i with respect to its equilibrium 
position, and Pi its momentum. The stiffness of the harmonic spring and the lattice 
spacing have been set to one, without losing generality. The coefficient a <C 1 measures 
the strength of the nonlinear contribution to the interaction potential. The two ends 
of the chain were assumed to be fixed, i.e. uq = ujv+i = 0. 

§ To whom correspondence should be addressed (Thierry.Dauxois@ens-lyon.fr) 

I Mary Tsingou, who took part in the numerical study is not an author of the report, but her 
contribution is however recognised by two lines of acknowledgments! 
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The common approach in physics is to think in terms of "normal modes" , 

related to the displacements through A}^ — ^2/ {N + 1) X]i3=i sin {ikir/N + 1) with 

the frequencies uj^ = 4sin^ (fc7r/2(A^ + 1)). They provide a basis to rewrite the 
Hamiltonian as 

^ N N 

-ff = 2 X! (^fc + ^fe^fe) + ^ X! CkemAkAiAmUJk(^euJm , (2) 

fc=l k,l,m=l 

where the coefficients Ckim are given for example in Ref. The last term, generated 
by the nonlinear contribution to the potential, leads to a coupling between the modes, 
and scales as N^^"^. 




Figure 1. FPU recurrence: the plot shows the time evolution of the kinetic 
and potential energy Ek = \ (A| + u/^A^ of each of the four lowest modes. 
Initially, only mode 1 was excited (from reference pQ). 

FPU thought that, due to this term, the energy introduced into a single mode, 
mode fc = I in their simulation, should slowly drift to the other modes, until the 
equipartition of energy predicted by statistical physics is reached. The beginning of 
the calculation indeed suggested that this would be the case. Modes 2,3,..., were 
successively excited. However, by accident 0], one day, they let the program run long 
after the steady state had been reached. When they realized their oversight and came 
back to the room, they noticed that the system, after remaining in a steady state 
for a while, had then departed from it. To their great surprise, after 157 periods of 
the mode k — 1, almost all the energy (all but 3%) was back to the lowest frequency 
mode, as shown in figure ^ 

The initial state seems to be almost perfectly recovered after this recurrence 
period. Further calculations, performed later with faster computers, showed that the 
same phenomenon repeats many times, and that a "super-recurrence" period exists, 
after which the initial state is recovered with a much higher accuracy. Thus, contrary 
to the expectations of the authors, the drift of the energy of the initially excited 
mode 1 towards the modes with a higher wave number does not occur. This highly 
remarkable result, known as the FPU paradox, shows that nonlinearity is not enough 
to guarantee equipartition of energy. 
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2. Fermi, Pasta, Ulam: the characters. 

Born in Rome, Enrico Fermi (1901-1954) has been one of the brightest physicists 
of the twentieth century, who made major experimental and theoretical discoveries. 
His name is famous for his contributions to statistical physics, elementary particle 
physics, the control of nuclear energy. It is thanks to his trip to receive the Nobel 
prize in Sweden in 1938 that he left fascist Italy, and emigrated to the United States 
where he studied atomic fission and set up the first controlled chain reaction in Chicago 
in 1942. He was naturally called to be one of the leaders of the Manhattan project for 
the development of nuclear energy and the atomic bomb. We are less familiar with his 
work on various nonlinear problems at the end of his life, before his premature death 
from stomach cancer. 

John R. Pasta (1918-1981) did not have such a bright career. In spite of 
his interest in physics he had to leave the New York City College during the great 
depression. He was first a real estate agent and then a police officer in New York 
from 1941 to 1942. Then he has been recruited as radar and cryptography specialist 
by the American army during the second world war. His earnings as a GI allowed 
him to return to the University, where he got a PhD in theoretical physics in 1951. 
He immediately started to work in Los Alamos, on the MANIAC computer which 
was under construction and tests. His career went developed as a computer expert 
for the Atomic Energy Commission, where he extensively developed the Mathematics 
and Computer Division. Then John R. Pasta became a physics professor and, in 1964, 
dean of the computer science department of the university of Illinois, where he focused 
on the use of computers to solve applied problems in physics and mathematics. 

Born in Poland, Stanislaw M. Ulam (1909-1984) quickly became interested in 
mathematics and obtained a PhD in 1933 under the supervision of Banach (1892-1945). 
Following a first invitation by Von Neumann to the prestigious Princeton Institute for 
Advanced Study in 1935, he visited the Unites States several times before leaving 
Poland definitely before the start of the second world war. He became an American 
citizen in 1943, and was invited by von Neumann himself to join the Los Alamos 
team to prepare the atomic bomb. Among other things he solved the problem of the 
initiation of fusion in the hydrogen bomb. In collaboration with N. C. Metropolis he 
invented the Monte-Carlo method to solve problems requiring a statistical sampling. 
He stayed in New Mexico until 1965 before his appointment as a mathematics professor 
at Colorado University. 

The FPU numerical experiment has been performed on the MANIAC 
(Mathematical Analyser, Numerical Integrator And Computer) built in 1952 for the 
Manhattan project, which has been used in the development of Mike, the first hydrogen 
bomb. Richard Feynman (1918-1987) and Nicholas C. Metropofis (1915-1999), 
exasperated by the slow and noisy mechanical calculators which were nevertheless 
necessary to the design of the bombs, promoted its construction. The name of the 
computer, MANIAC, had been chosen by the project director, N. Metropolis, who 
hoped that it would stop the raising fashion of naming computers by acronyms. The 
effect has been exactly the opposite, although, according to Metropolis "5, George 
Gamow (1904-1968) was instrumental in rendering this and other computer names 
ridiculous when he dubbed the Maniac "Metropolis And von Neumann Install Awful 
Computer" . The Maniac was able to perform about 10** operations per second, which 
must be compared to the 10^ operations per second of any personal computer today. 
It is Fermi who had the genius to propose that computers could be used to study a 



The Fermi-Pasta- Ulam numerical experiment 



4 



problem or test a physical idea by simulation, instead of simply performing standard 
calculus. He proposed to check the prediction of statistical physics on this system now 
called FPU. 

The discovery that resulted from this study has not only been at the origin of 
the soliton concept and of many features of chaotic phenomena, as discussed in the 
following section, but it also introduced the concept of numerical experiment. This 
had far reaching consequences, leading to a complete revolution in the investigation of 
physical phenomena. The computer is no more used only to perform a calculation that 
cannot be done by pencil-and-napkin, but to check a theoretical conjecture that cannot 
be proven analytically, or even to provide the theorist with "experimental" results 
that wait for a mathematical proof: a source of problems, like in a "true laboratory 
experiment" . Of course, the numerical experiment does not possess all the complexity 
of a true physical experiment, because the reality which is represented is highly virtual. 
However, today numerical simulations of condensed matter systems, one reaches such 
a level of confidence that sometimes it has happened that a laboratory experiment has 
been questioned because a numerical experiment had given contrary indications 
Today, computational physics is an established discipline and it is considered as 
sort of separated from both theoretical and experimental physics. Students are 
currently trained in computational physics as in other disciplines, and specialized 
journals publish the results of the research in this field. This big epistemological and 
sociological change began with the FPU experiment. 

3. Epilogue: The solution of the FPU problem 

Pursuing the solution of the FPU paradox, two main different lines of thought were 
followed. On one side, people like J. Ford [Hj focused on the Fourier mode dynamics, 
looking for non-resonance conditions that could explain the inefficient energy transfer. 
No convincing explanation was found before the news spread of a theorem announced 
by Kolmogorov [7|, and then proven by Moser and Arnold |9| (KAM theorem), 
which states that most orbits of a slightly perturbed integrable Hamiltonian systems 
remain quasi-periodic. Although qualitative explanations of the FPU recurrence 
were obtained following this path no quantitative explanation has been obtained, 
to our knowledge. On the contrary, a remarkable result was obtained by Izrailev 
and Chirikov ^U], following KAM theory. If the perturbation is strong enough that 
nonlinear resonances "superpose" , the FPU recurrence is destroyed and one obtains a 
fast convergence to thermal equilibrium. This prediction was tested numerically and is 
nowadays at the basis of the phenomenon called "strong stochasticity threshold" ^Tj . 

The other line of thought, which went towards the so-called continuum limit, led 
to the solution of the FPU paradox by Zabusky and Kruskal in terms of the 
dynamics of solitons ^S]- Starting from the FPU equations of motion, derived from 
Hamiltonian QJ, 

ili = {ui+i -^Ui^i -2ui) a[(ui+i- - {ui- Ui-iY] , (3) 

and restricting the investigation to long wavelength modes, they could derive the 
Korteweg-de Vries equation ^1] 

Wr + + " = , (4) 

where the field w is a conveniently rescaled spatial derivative of the displacement field 
Ui, r is a rescaled time and Wx = dw/dx. 
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The sinusoidal initial condition develops sharp fronts, and then breaks into a 
series of pulses, which are solitons. They preserve their shapes and velocities and, 
during their motion in the finite system with periodic boundary conditions^, from 
time to time, they come back to the positions that they had initially, restoring the 
initial condition. This is why FPU recurrence is observed. 

It is important to emphasize that the FPU paradox would probably not have been 
a mystery for more than 10 years, if, before Zabusky and Kruskal, somebody had the 
idea to look carefully at the dynamics of the nonlinear lattice as a function of the space 
coordinate (see the suggested numerical experiment in Sec. ISjl. But physicists were 
so used to analyse linearised problems, which naturally leads to a description in the 
Fourier space, and only add nonlinearity afterwards as a coupling between the modes, 
that the observation of the solitons emerging from the sinusoidal initial condition had 
not been done. 

4. The remake: FPU and the Japanese School 

The history of the FPU problem is quite rich, as documented in the report by Ford ^H] 
and in the recent book by Weissert ^^I- We will not discuss it in this paper and we 
just recall that some of the contributions to the Chaos issue JZj, printed this year 
to celebrate the anniversary of the FPU experiment, are devoted to a brief historical 
reconstruction. 

Because of its relevance for the determination of a quantitative formula for the 
recurrence time, we would like to discuss briefly the contribution of the Japanese 
school, which has been very active in this field since the 60's. 

The original FPU report (J, not published, was known by very few people, most 
of them within America. This is why, in Japan, Nobuhiko Saito (1919-) together 
with his PhD student Hajime Hirooka (1939-) were developing closely related research 
works in complete ignorance of what had been done in Los Alamos ten years earlier. 

As explained by the title of Hirooka's PhD Thesis, "The approach to thermal 
equilibrium in a nonlinear lattice" , the motivation was to explore the mechanism of 
ergodicity. As the mechanics of collisions in a gas was too complicated, they thought 
that an anharmonic lattice would be a more appropriate system; it appeared finally 
that its dynamics was difficult to solve. This is why, in 1964 they started to perform 
numerical simulations on a small computer provided by NEC: it took all the night 
to compute the evolution for only 5 lattice sites during several periods of the lowest 
mode! In 1965, they switched to the new Supercomputer provided by Hitachi Co. 
to Tokyo University, installed under a national plan for all Japanese universities and 
research institutes. 

Saito and Hirooka considered a one-dimensional anharmonic lattice with 
quadratic and quartic potentials between neighboring particles, and with both ends 
fixed. The initial excitation was slightly different since all particles were at rest, 
whereas a constant force was applied to the first particle. They also prepared a similar 
system, with only harmonic potentials, and analytically calculated several quantities 
of interest, in particular the long-time averages of the squares of the velocities of the 
particles {u^) and the correlations {iiiiii+i) of the velocities of neighboring particles. 
The long-time averages of (uf) were the same for all particles, but the correlation 

^ The calculation performed by Zabusky and Kruskal used periodic boundary conditions, while the 
FPU calculation used fixed boundaries. This does not change the analysis because either the solitons 
pass through the boundaries and enter on the opposite side, or they are totally reflected 
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functions did not vanish in the long term, and the MaxweUian distribution of velocities 
was not observed, contrary to their expectations. 

It is interesting to emphasize that they were hesitating to publish these results, 
because they were not familiar with the computer calculation and were afraid of 
having introduced some errors. However, they found numerical experiments similar 
to the FPU original one in several papers on nonlinear oscillation, in particular, 
those of Ford [Hj. These papers also convinced Morikazu Toda to study solitons 
and to introduce the lattice with exponential interactions, nowadays called the Toda 
lattice [IHl. 

Soon after, they also learned that the original FPU report was reproduced 
in the Collected Papers of Fermi. Very excited, Saito and Toda read the paper at 
the library of Tokyo University of Education (now Tsukuba University). As they had 
found results similar to those of Fermi, Pasta and Ulam, they finally decided to publish 
their calculations ^^l- Afterwards, they considered the simpler FPU initial condition 
and found the induction phenomenon and the occurrence of the random character of 
lattice vibrations. 

Finally, they discovered the Zabusky-Kruskal seminal paper soon after the 
publication of their papers in 1967. Norman Zabusky went to Kyoto in 1968 at 
the International conference of Statistical Mechanics, where he showed his movies 
and in particular the KdV cinema, which "has the power to communicate unbiased 
information in a credible manner, much beyond the power of words, graphs and 
equations" . 

At that time, M. Toda gave the first analytical estimate of the recurrence time, 
based on the exact solutions of the exponential lattice, nowadays called the Toda 
lattice. Introducing the period Ti = 2N of the first mode and the amplitude a of the 
initial sine wave, he got 

3 iV3/2 ^5/2 

Tr = ^ ^ Ti ~ 0.76 ^= , (5) 



which compares well with the first empirical estimate due to Zabusky. Toda's result 
was a little larger than Zabusky's formula: the prefactor is 0.31 instead of 0.38 in 
equation ((SJ. The discrepancy is due to the fact that, when solitons pass through each 
other, they accelerate because of the compression of the lattice. Sec Ref. j20j for a 
more complete discussion and Ref. 21^ for a recent check of formula l(SJl. 



5. Pedagogical perspectives 

It is quite easy, with contemporary computational tools, to reproduce the original FPU 
experiment. An example is the simple MATLAB© code reported in the Appendix, 
which solves the equations for the dynamics of the FPU model and computes the 
normal modes Ak- We suggest the student to use it to reproduce different aspects of 
the FPU problem 

• Try to reproduce Fig. fusing the original FPU initial condition given in the code. 
Be careful in the choice of the amplitude a because the cubic lattice is unstable 
at large energies. 

• By careful inspection of the difference field (u^+i — Ui), it is possible to detect the 
presence of propagating structures, which are indeed the solitons of the Zabusky- 
Kruskal analysis. 
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• Between a — \ and a — 10, the FPU model with = 32 undergoes relaxation to 
equilibrium ^U] . Check the absence of FPU recurrences and the establishment of 
energy equipartition among the normal modes. 

• Since the initial sine wave form solitons |12| . it is quite natural to begin with a 
lattice soliton initial condition. The exact kink-like soliton solution for the Toda 
lattice EDI is 



where k is the inverse width. This solution is only approximate for the FPU 
lattice, and due to the fixed boundary conditions, one must put a kink (+) and 
an antikink (-) solutions, centered around 2 different sites of the lattice. Our 
advice is to first choose a large lattice (e.g. N — 128). The time evolution of 
the derivative field Ui+i — Ui shows the propagation of coherent structures, their 
mutual interactions and reflection from the boundaries. 

• Try to reproduce the relevant scalings of formula ((S)), for instance the N'^^^ 
dependence on the number of oscillators. 

• Finally, try the Zabusky-Deem [221 initial condition Ui — asin[i7r7V/(iV + 1)], 
which corresponds to the highest frequency mode. Similar oscillations of mode 
energies as the one of long wavelengths are observed. Try to characterize these 
recurrences. 

The scaling of formula Q is quite different from the one of the Poincare recurrence 
time |2,S| . which can be shown to increase exponentially with N in the harmonic 
limit 123]. No striking change of this scaling should occur for a small nonlinearity, 
e.g. in the FPU recurrence region, although nothing is rigorously known for this 
case. The interested student can find the statement of Poincare recurrence theorem 
in most books of statistical mechanics. However, in simple words, the theorem says 
that whenever a dynamical system preserves phase-space volume in its time evolution 
(and this is the case of Hamiltonian systems, including the FPU oscillator chain), 
"almost all" trajectories return arbitrarily close to their initial position, and they do 
this an infinite number of times. The exponential scaling with N of the Poincare 
recurrence time with respect to the algebraic one of Toda's formula proves that the 
FPU recurrence is not a manifestation of this phenomenon. It's remarkable that, 
although certainly Fermi and Ulam did know Poincare theorem, they did not quote it 
at all in their paper as a possible explanation of the observed recurrence"*". 

Another pedagogical aspect, which is worth discussing in connection with the 
FPU experiment, is the relation between continuum and discrete. Analytical "human" 
approaches to the problem consider continuous variables, while with the computers, 
it was natural and necessary to use discretized variables; the human and computer 
way of dealing with physical problems were therefore for the first time contrasted. 
The only discreteness appearing in the FPU model is the spatial one: the oscillators 
are put on a lattice. Time flows continuously and the displacement field variable 
Ui is also continuous. Models of space-time evolution have been considered where 
time is discrete, while the lattice variable is still continuous: so-called coupled map 
lattices [21]. Hamiltonian versions of this type of models exist, in the sense that phase- 
space volume is conserved by the dynamical evolution, although the total energy is not. 

+ This fallacious possibility has been considered by others (who perhaps did not know the theorem 
as well as Fermi and Ulam) when the news about the FPU experiment spread around. 



Ui{t) = ±— log 



1 -I- exp[2fc(i — 1 — io) ± i sinh k] 



(6) 



1 -I- exp[2A:(i — io) ± i sinh k] 
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The extreme case of discreteness is the one of cellular automata |26| . for which also 
the lattice variable takes a discrete set of values. Depending on the physical system 
at hand, one model or the other has been considered as appropriate, and several 
examples exist where even the extreme discrete case of cellular automata provides a 
good representation of the physical process under analysis. 

In this respect, quite advanced still open questions related to the FPU recurrence 
are: Is the recurrence observed also in "Hamiltonian" coupled map lattices and in 
cellular automata? Is it related to the propagation of coherent structures which have 
some resemblance to solitons? How does the recurrence time scale with system size? 
These are open research issues, on which perhaps a curious PhD student might want 
to challenge his own computational skill, taking the opportunity to make his own first 
"numerical experiment" , a privilege once reserved only to a few people isolated in a 
deserted region of New Mexico. 
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Appendix: MATLAB© Code 

Codes are available at the website: 

[http : //perso . ens-lyon . f r/Thierry . Dauxois/f pu . htmlj 

Main Program 

N=32; °/„ Number of particles must be a power of 2 

alpha=0.25; 7, Nonlinear parameter 

TMAX=10000; DT=20; 7. tmax and Delta t 

tspan=[0:DT:TMAX] ; 

options=odeset ('Reltol' , le-4, 'DutputFcn' , 'odeplot ' , 'OutputSel' , [1,2,N] ) ; 
7o Test different tolerances, changing Reltol 
for 1=1:N, 

a=l; b(l)=a*sin(pi*l/(N+l)) ; b(l+N)=0; 7. FPU initial condition 

7.a=l; b(l)=a*sin(pi*N*l/(N+l)) ; b(l+N)=0; 7. Zabusky-Deem init . cond. 
7.k=0.8; sk=(sinh(k))"2; ek=exp(-k); il=l-N/4; i2=il-N/2; 7.Solitons 
7.b(l)=-0.5/alpha*log((l+exp(2*k*(il-l)))/(l+exp(2*k*il))) ; 
7.b(l)=b(I)+0 . 5/alpha*log( (l+exp(2*k* (i2-l) ) ) / (l+exp(2*k*i2) ) ) ; 
7.b(I+N)= sk*ek/alpha/cosh(k*il)/(exp(-k*il)+exp(k*il)*exp(-2*k)) ; 
7.b (1+N) =b (1+N) -sk*ek/alpha/cosh (k*i2) / (exp (-k*i2) +exp (k*i2) *exp (-2*k) ) ; 
omegak2(l)=4*(sin(pi*l/2/N))"2; 7. Mode Frequencies 



end 



[T,Y]=ode45('fpul' ,tspan,b' , options, N) ; 
for 1T=1: (TMAX/DT) , 



7o Time integration 
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TIME(IT)=IT*DT*sqrt(oinegak2(l))/2/pi; 7, Time iteration loop 

YX(IT,1:N+1) = [0 Y(IT,1:N )] ; YV(IT, 1 : N+l) = [0 Y(IT,N+1:2*N )] ; 

sXF ( IT , : ) =imag (f f t ( [YX ( IT , 1 : N+1) -YX (IT , N+1 : - 1 : 2) ] ) ) /sqrt (2* (N+1 ) ) ; 

sVF(IT, :)=imag(fft([YV(IT,l:N+l) -YV(IT,N+1 : -1 : 2)] ) )/sqrt (2* (N+1) ) ; 

Energ(IT,l:N)=(omegak2(l:N) .*(sXF(IT,2:N+l) ."2)+sVF(IT,2:N+l) .~2)/2; 

for J=2:N-1, 7. Space loop 

DifYdT, J)=Y(IT, J+1)-Y(IT, J) ; 

end 

end 

plot(TIME,Energ(: ,1) ,TIME,Energ( : ,2) ,TIME,Energ( : ,3) ,TIME,Energ(: ,4)) ; 
surf(DifY); '/„ Space derivative field to show the soliton dynamics 

fpul function 

function dy=f pul (t ,y) ; 
N=32 ; alplia=0 . 25 ; 

D(N+l)=y(2) -2*y(l)+alpha*((y(2)-y(l))-2-y(l)-2) ;D(l)=y(N+l) ; 
D(2*N)=y(N-l)-2*y(N)+alpha*(y(N)'2-(y(N)-y(N-l))~2) ;D(N)=y(2*N) ; 
for I=2:N-1, 

D(N+I)=y(I+l)+y(I-l)-2*y(I)+alpha*((y(I+l)-y(I))-2-(y(I)-y(I-l))-2); 

D(I)=y(N+I) ; 

end 

dy=D'; 
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